First, check package dependencies and install ProjecTILs
Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true")
remotes::install_github("carmonalab/UCell", ref="v1.3")
remotes::install_github("carmonalab/scGate")
remotes::install_github("carmonalab/ProjecTILs")
library(ProjecTILs)
library(Seurat)
First, load the default reference TIL atlas. If no reference map file
is provided, the function load.reference.map() will
automatically download it from https://doi.org/10.6084/m9.figshare.12478571
ref <- load.reference.map()
## [1] "Loading Default Reference Atlas..."
## [1] "/Users/mass/Documents/Projects/Cell_clustering/ProjecTILs.demo/ref_TILAtlas_mouse_v1.rds"
## [1] "Loaded Reference map ref_TILAtlas_mouse_v1"
Let’s explore the reference atlas
refCols <- c("#edbe2a", "#A58AFF", "#53B400", "#F8766D", "#00B6EB", "#d1cfcc", "#FF0000", "#87f6a5", "#e812dd")
DimPlot(ref,label = T, cols = refCols)
See expression of important marker genes across reference subtypes
markers <- c("Cd4","Cd8a","Ccr7","Tcf7","Pdcd1","Havcr2","Tox","Izumo1r","Cxcr6","Xcl1","Gzmb","Gzmk","Ifng","Foxp3")
VlnPlot(ref,features=markers,stack = T,flip = T,assay = "RNA")
Now let’s load a query dataset - Miller et al., Nature Immunol (2019)
#A sample data set is provided with the ProjecTILs package
querydata <- ProjecTILs::query_example_seurat
More generally, it is possible to load a query matrix with gene names and barcodes (e.g. 10X format or raw counts)
##Raw count matrix from GEO
library(GEOquery)
geo_acc <- "GSE86028"
getGEOSuppFiles(geo_acc)
fname3 <- sprintf("%s/GSE86028_TILs_sc_wt_mtko.tpm.log2.txt.gz", geo_acc)
querydata3 <- read.sc.query(fname3, type = "raw.log2")
query.projected <- make.projection(querydata, ref=ref)
## [1] "Using assay RNA for query"
## [1] "37 out of 1501 ( 2% ) non-pure cells removed. Use filter.cells=FALSE to avoid pre-filtering"
## [1] "Aligning query to reference map for batch-correction..."
##
## Projecting corrected query onto Reference PCA space
##
## Projecting corrected query onto Reference UMAP space
NB: by default, make.projection() will pre-filter T
cells using scGate.
In case the input dataset is already pre-filtered, or if you are using a
non-T cell reference atlas, you can disable this step using
make.projection(querydata, ref=ref, filter.cells = FALSE).
Plot projection of new data over the reference in UMAP space. The contour lines display the density of projected query cells onto the reference map.
plot.projection(ref, query.projected)
Predict the cell states in the query set using a nearest-neighbor algorithm
query.projected <- cellstate.predict(ref=ref, query=query.projected)
table(query.projected$functional.cluster)
##
## CD8_EffectorMemory CD8_NaiveLike CD8_Tex CD8_Tpex
## 61 24 1308 58
## Th1 Treg
## 1 12
Plot the predicted composition of the query in terms of reference T cell subtypes
plot.statepred.composition(ref, query.projected,metric = "Percent")
How do the gene expression levels compare between reference and query for the different cell states?
plot.states.radar(ref, query=query.projected, min.cells=30)
If we have multiple conditions (e.g. control vs. treatment, or samples from different tissues), we can search for discriminant genes between conditions (otherwise, by default this analysis is performed against the reference subtype as the ‘control’)
#Simulate a condition which e.g. increases Gzmb expression compared to control
query.control <- subset(query.projected, subset=`Gzmb` < 1.5)
query.perturb <- subset(query.projected, subset=`Gzmb` >= 1.5)
plot.states.radar(ref, query=list("Control" = query.control, "Query" = query.perturb))
In this toy example, where we simulated a condition that increases Gzmb expression compared to control, we expect cytotoxicity genes to drive differences.
discriminantGenes <- find.discriminant.genes(ref=ref, query=query.perturb, query.control=query.control, state="CD8_Tex")
head(discriminantGenes,n=10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Gzmb 3.842405e-94 2.8257408 1.000 0.848 1.050552e-89
## Gzmc 2.615333e-35 2.8814254 0.809 0.378 7.150581e-31
## Lgals1 5.268875e-24 0.4543825 1.000 1.000 1.440563e-19
## AA467197 2.264136e-22 1.1983647 0.845 0.530 6.190373e-18
## Mt1 3.582293e-19 1.1840434 0.776 0.457 9.794346e-15
## Ccl9 1.921541e-18 1.6471873 0.489 0.116 5.253686e-14
## Serpinb6b 3.538035e-16 0.8050444 0.858 0.598 9.673341e-12
## Mt2 1.750597e-14 0.9658675 0.613 0.274 4.786308e-10
## Gzme 4.430793e-14 2.0820514 0.462 0.165 1.211423e-09
## Prf1 5.926832e-14 0.5883041 0.503 0.177 1.620455e-09
We can use a volcano plot to display differentially expressed genes:
library(EnhancedVolcano)
## Loading required package: ggrepel
## Registered S3 methods overwritten by 'ggalt':
## method from
## grid.draw.absoluteGrob ggplot2
## grobHeight.absoluteGrob ggplot2
## grobWidth.absoluteGrob ggplot2
## grobX.absoluteGrob ggplot2
## grobY.absoluteGrob ggplot2
EnhancedVolcano(discriminantGenes, lab = rownames(discriminantGenes), x = "avg_log2FC", y = "p_val", pCutoff = 1e-09,
FCcutoff = 0.5, labSize = 5, legendPosition = "none", drawConnectors = F, title = "Gzmb_high vs. Gzmb_low (Tex)")
## Warning: Ignoring unknown parameters: xlim, ylim
Using a random subsetting, p-values should not be significant:
rand.list <- ProjecTILs:::randomSplit(query.projected, n=2, seed=1)
discriminantGenes <- find.discriminant.genes(ref=ref, query=rand.list[[1]], query.control=rand.list[[2]], state="CD8_Tex")
EnhancedVolcano(discriminantGenes, lab = rownames(discriminantGenes), x = "avg_log2FC", y = "p_val", pCutoff = 1e-09,
FCcutoff = 0.5, labSize = 5, legendPosition = "none", drawConnectors = F, title = "Random split (Tex)")
The dimensions in UMAP space summarize the main axes of variability of the reference map. What if the query data contains novel states? We can search for additional, maximally discriminant dimensions (either in ICA or PCA space) that explain new variability in the query set.
As before, simulate a condition which increases Gzmb expression compared to control
#
query.control <- subset(query.projected, subset=`Gzmb` < 1.5)
query.perturb <- subset(query.projected, subset=`Gzmb` >= 1.5)
plot.states.radar(ref, query=list("Control" = query.control, "Query" = query.perturb))
In this toy example, we expect some gene module associated with granzymes to drive the discriminant analysis:
top.ica.wcontrol <- find.discriminant.dimensions(ref=ref, query=query.perturb, query.control=query.control)
head(top.ica.wcontrol)
## stat stat_abs p_val
## ICA_26 0.4800657 0.4800657 0.000000e+00
## ICA_24 -0.3245565 0.3245565 7.521761e-12
## ICA_19 0.2854554 0.2854554 7.046419e-09
## ICA_42 0.2608306 0.2608306 3.341548e-07
## ICA_6 0.2407044 0.2407044 6.047000e-06
## ICA_28 -0.2129243 0.2129243 2.246895e-04
VizDimLoadings(ref, reduction = "ica", nfeatures = 10, dims=c(26,24,42), ncol=3)
Now we can plot the ICA dimension that captured the genetic changes associated to the perturbation of increasing Gzmb
plot3d <- plot.discriminant.3d(ref, query=query.perturb, query.control=query.control, extra.dim="ICA_26")
plot3d
We can plot other metadata in the z-axis of the UMAP, e.g. the cycling score calculated by the TILPRED cycling signature
plot3d <- plot.discriminant.3d(ref, query.projected, extra.dim="cycling.score")
plot3d
Focus the plot only on a specific state
plot3d <- plot.discriminant.3d(ref, query.projected, extra.dim="cycling.score", query.state="CD8_Tex")
plot3d
Using a random subsetting, p-values should not be significant:
rand.list <- ProjecTILs:::randomSplit(query.projected, n=2, seed=1)
top.ica.ks.rand <- find.discriminant.dimensions(ref=ref, query=rand.list[[1]], query.control=rand.list[[2]], reduction="ica")
top.ica.ttest.rand <- find.discriminant.dimensions(ref=ref, query=rand.list[[1]], query.control=rand.list[[2]], reduction="ica", test = "t-test")
ProjecTILs repository
ProjecTILs case studies - INDEX - Repository
Publication: Andreatta et al Nat. Comm. 2021